"""
DaoMath Thermodynamics & Statistical Physics (道数热力学)
从 ±0 先天均衡构造热力学与统计力学

热力学的 ±0 诠释:
    吸热 ↔ 放热 (对偶)
    有序 ↔ 无序 (对偶)
    熵增 ↔ 能量守恒 (对偶)
    
理论层级:
    ±0 (先天均衡)
     ↓
    微观态 (Microstates)
     ↓
    统计系综 (Statistical Ensemble)
     ↓
    配分函数 (Partition Function)
     ↓
    热力学量 (Thermodynamic Quantities)
     ↓
    热力学定律
     ↓
    相变与临界现象
     ↓
    非平衡态
"""

from __future__ import annotations
from typing import (
    List, Tuple, Callable, Optional, Dict, Any,
    Union, TypeVar, Generic
)
from dataclasses import dataclass, field
from abc import ABC, abstractmethod
from enum import Enum
import numpy as np
import scipy.linalg as la
from scipy.integrate import quad, odeint
from scipy.optimize import minimize, fsolve
from scipy.special import factorial, comb
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

from core.axioms import InnateBalance, DualPair


# ============================================================================
# 物理常数
# ============================================================================

class PhysicalConstants:
    """物理常数"""
    K_B = 1.380649e-23          # Boltzmann 常数 (J/K)
    N_A = 6.02214076e23         # Avogadro 常数 (1/mol)
    R = 8.314462618             # 气体常数 (J/(mol·K))
    H = 6.62607015e-34          # Planck 常数 (J·s)
    HBAR = 1.054571817e-34      # 约化 Planck 常数
    C = 299792458.0             # 光速 (m/s)
    SIGMA = 5.670374419e-8      # Stefan-Boltzmann 常数 (W/(m²·K⁴))
    
    # 简化单位制
    K_B_REDUCED = 1.0


# ============================================================================
# Level 1: 微观态与统计系综
# ============================================================================

class Microstate:
    """
    微观态
    
    系统的一个具体量子态或相空间点
    
    基于 ±0 的理解:
    - 每个微观态：一种可能性
    - 宏观态：微观态的集合（正负均衡）
    """
    
    def __init__(self, 
                 energy: float,
                 quantum_numbers: Optional[List[int]] = None,
                 degeneracy: int = 1):
        """
        Args:
            energy: 能量本征值
            quantum_numbers: 量子数 (n, l, m, s, ...)
            degeneracy: 简并度
        """
        self.energy = energy
        self.quantum_numbers = quantum_numbers or []
        self.degeneracy = degeneracy
    
    def __repr__(self) -> str:
        return f"Microstate(E={self.energy:.4f}, g={self.degeneracy})"


class EnergyLevel:
    """
    能级
    
    具有相同能量的微观态集合
    """
    
    def __init__(self, energy: float, degeneracy: int):
        self.energy = energy
        self.degeneracy = degeneracy
    
    def boltzmann_factor(self, temperature: float, 
                        k_B: float = PhysicalConstants.K_B_REDUCED) -> float:
        """
        Boltzmann 因子: e^(-E/k_B T)
        """
        beta = 1.0 / (k_B * temperature)
        return np.exp(-beta * self.energy)


class StatisticalEnsemble(ABC):
    """
    统计系综抽象基类
    
    系综：大量相同宏观条件的系统集合
    
    类型:
    - 微正则系综: 孤立系统 (E, V, N 固定)
    - 正则系综: 恒温系统 (T, V, N 固定)
    - 巨正则系综: 开放系统 (T, V, μ 固定)
    """
    
    @abstractmethod
    def partition_function(self) -> float:
        """配分函数"""
        pass
    
    @abstractmethod
    def average_energy(self) -> float:
        """平均能量"""
        pass
    
    @abstractmethod
    def entropy(self) -> float:
        """熵"""
        pass


class MicrocanonicalEnsemble(StatisticalEnsemble):
    """
    微正则系综
    
    孤立系统: E, V, N 固定
    
    基本假设: 所有可达微观态等概率
    """
    
    def __init__(self, 
                 energy: float,
                 volume: float,
                 particle_number: int,
                 microstates: List[Microstate]):
        """
        Args:
            energy: 总能量
            volume: 体积
            particle_number: 粒子数
            microstates: 满足条件的微观态列表
        """
        self.E = energy
        self.V = volume
        self.N = particle_number
        self.microstates = microstates
        self.Omega = len(microstates)  # 微观态数目
    
    def partition_function(self) -> float:
        """
        Ω(E, V, N): 微观态数
        
        这不是传统意义的配分函数，而是状态数
        """
        return float(self.Omega)
    
    def average_energy(self) -> float:
        """平均能量（就是给定的能量）"""
        return self.E
    
    def entropy(self, k_B: float = PhysicalConstants.K_B_REDUCED) -> float:
        """
        Boltzmann 熵公式
        
        S = k_B ln Ω
        
        体现 ±0: 最大无序（均衡）对应最大熵
        """
        if self.Omega <= 0:
            return 0.0
        return k_B * np.log(self.Omega)
    
    def temperature(self, k_B: float = PhysicalConstants.K_B_REDUCED) -> float:
        """
        温度定义
        
        1/T = (∂S/∂E)_{V,N}
        """
        # 数值微分
        dE = 1e-6
        
        # 计算 E + dE 时的熵（需要新的微观态数）
        # 简化：假设 Ω ∝ E^(3N/2) (理想气体)
        dOmega_dE = (3 * self.N / 2) * self.Omega / self.E if self.E > 0 else 0
        
        dS_dE = k_B * dOmega_dE / self.Omega if self.Omega > 0 else 0
        
        if dS_dE <= 0:
            return float('inf')
        
        return 1.0 / dS_dE


class CanonicalEnsemble(StatisticalEnsemble):
    """
    正则系综
    
    恒温系统: T, V, N 固定
    
    与热库接触的系统
    """
    
    def __init__(self,
                 temperature: float,
                 volume: float,
                 particle_number: int,
                 energy_levels: List[EnergyLevel]):
        """
        Args:
            temperature: 温度 T
            volume: 体积 V
            particle_number: 粒子数 N
            energy_levels: 能级列表
        """
        self.T = temperature
        self.V = volume
        self.N = particle_number
        self.energy_levels = energy_levels
        self.k_B = PhysicalConstants.K_B_REDUCED
        self.beta = 1.0 / (self.k_B * self.T)
    
    def partition_function(self) -> float:
        """
        正则配分函数
        
        Z = Σ_i g_i e^(-β E_i)
        
        体现 ±0: 低能态（负指数大）vs 高能态（正指数小）的均衡
        """
        Z = 0.0
        for level in self.energy_levels:
            Z += level.degeneracy * np.exp(-self.beta * level.energy)
        return Z
    
    def average_energy(self) -> float:
        """
        平均能量
        
        ⟨E⟩ = -∂ln Z/∂β = Σ_i E_i P_i
        """
        Z = self.partition_function()
        if Z == 0:
            return 0.0
        
        E_avg = 0.0
        for level in self.energy_levels:
            prob = level.degeneracy * np.exp(-self.beta * level.energy) / Z
            E_avg += level.energy * prob
        
        return E_avg
    
    def entropy(self) -> float:
        """
        熵
        
        S = k_B (ln Z + β⟨E⟩)
        """
        Z = self.partition_function()
        if Z <= 0:
            return 0.0
        
        E_avg = self.average_energy()
        return self.k_B * (np.log(Z) + self.beta * E_avg)
    
    def free_energy(self) -> float:
        """
        Helmholtz 自由能
        
        F = -k_B T ln Z = ⟨E⟩ - TS
        
        体现能量-熵的对偶
        """
        Z = self.partition_function()
        if Z <= 0:
            return float('inf')
        
        return -self.k_B * self.T * np.log(Z)
    
    def heat_capacity(self) -> float:
        """
        热容
        
        C_V = ∂⟨E⟩/∂T = k_B β² ⟨(E - ⟨E⟩)²⟩
        """
        Z = self.partition_function()
        if Z == 0:
            return 0.0
        
        E_avg = self.average_energy()
        
        # 计算 ⟨E²⟩
        E2_avg = 0.0
        for level in self.energy_levels:
            prob = level.degeneracy * np.exp(-self.beta * level.energy) / Z
            E2_avg += level.energy**2 * prob
        
        # 方差
        variance = E2_avg - E_avg**2
        
        return self.k_B * self.beta**2 * variance
    
    def probability_distribution(self) -> np.ndarray:
        """
        能级占据概率分布
        
        P_i = g_i e^(-β E_i) / Z
        """
        Z = self.partition_function()
        if Z == 0:
            return np.zeros(len(self.energy_levels))
        
        probs = np.array([
            level.degeneracy * np.exp(-self.beta * level.energy) / Z
            for level in self.energy_levels
        ])
        
        return probs


class GrandCanonicalEnsemble(StatisticalEnsemble):
    """
    巨正则系综
    
    开放系统: T, V, μ 固定
    
    与粒子库和热库接触
    """
    
    def __init__(self,
                 temperature: float,
                 volume: float,
                 chemical_potential: float,
                 single_particle_states: List[EnergyLevel]):
        """
        Args:
            temperature: 温度 T
            volume: 体积 V
            chemical_potential: 化学势 μ
            single_particle_states: 单粒子态
        """
        self.T = temperature
        self.V = volume
        self.mu = chemical_potential
        self.states = single_particle_states
        self.k_B = PhysicalConstants.K_B_REDUCED
        self.beta = 1.0 / (self.k_B * self.T)
    
    def partition_function(self) -> float:
        """
        巨正则配分函数
        
        Ξ = Σ_N Z_N e^(β μ N)
        """
        # 简化：对于非相互作用粒子
        Xi = 1.0
        for state in self.states:
            # 每个态可以有 0, 1, 2, ... 个粒子
            # 对于玻色子: Π_i 1/(1 - e^(-β(ε_i - μ)))
            # 对于费米子: Π_i (1 + e^(-β(ε_i - μ)))
            
            # 假设玻色子
            factor = -self.beta * (state.energy - self.mu)
            if factor < 0:
                Xi *= 1.0 / (1.0 - np.exp(factor))
            else:
                Xi *= 1.0  # 避免发散
        
        return Xi
    
    def average_particle_number(self) -> float:
        """
        平均粒子数
        
        ⟨N⟩ = k_B T ∂ln Ξ/∂μ
        """
        # 数值微分
        dmu = 1e-6
        Xi = self.partition_function()
        
        self.mu += dmu
        Xi_plus = self.partition_function()
        self.mu -= dmu
        
        if Xi <= 0 or Xi_plus <= 0:
            return 0.0
        
        dln_Xi_dmu = (np.log(Xi_plus) - np.log(Xi)) / dmu
        
        return self.k_B * self.T * dln_Xi_dmu
    
    def average_energy(self) -> float:
        """平均能量"""
        # ⟨E⟩ = -∂ln Ξ/∂β + μ⟨N⟩
        Xi = self.partition_function()
        
        # 简化计算
        E_avg = 0.0
        for state in self.states:
            # 玻色-Einstein 分布或费米-Dirac 分布
            occupation = self._occupation_number(state.energy)
            E_avg += state.energy * occupation
        
        return E_avg
    
    def entropy(self) -> float:
        """熵"""
        Xi = self.partition_function()
        if Xi <= 0:
            return 0.0
        
        E_avg = self.average_energy()
        N_avg = self.average_particle_number()
        
        return self.k_B * (np.log(Xi) + self.beta * (E_avg - self.mu * N_avg))
    
    def _occupation_number(self, energy: float) -> float:
        """单态占据数（玻色子）"""
        factor = self.beta * (energy - self.mu)
        if factor > 100:  # 避免数值溢出
            return 0.0
        
        denominator = np.exp(factor) - 1
        if abs(denominator) < 1e-10:
            return float('inf')
        
        return 1.0 / denominator


# ============================================================================
# Level 2: 理想气体模型
# ============================================================================

class IdealGas:
    """
    理想气体
    
    最简单的热力学模型
    
    基于 ±0: 粒子间无相互作用（净作用为零）
    """
    
    def __init__(self, 
                 n_particles: int,
                 temperature: float,
                 volume: float,
                 mass: float = 1.0):
        """
        Args:
            n_particles: 粒子数 N
            temperature: 温度 T
            volume: 体积 V
            mass: 粒子质量 m
        """
        self.N = n_particles
        self.T = temperature
        self.V = volume
        self.m = mass
        self.k_B = PhysicalConstants.K_B_REDUCED
    
    def pressure(self) -> float:
        """
        理想气体状态方程
        
        PV = Nk_B T
        """
        return self.N * self.k_B * self.T / self.V
    
    def internal_energy(self) -> float:
        """
        内能（单原子理想气体）
        
        U = (3/2) Nk_B T
        """
        return 1.5 * self.N * self.k_B * self.T
    
    def entropy(self) -> float:
        """
        Sackur-Tetrode 方程
        
        S = Nk_B [ln(V/N(mk_BT/2πℏ²)^(3/2)) + 5/2]
        """
        h = PhysicalConstants.H
        
        # 热波长
        lambda_th = h / np.sqrt(2 * np.pi * self.m * self.k_B * self.T)
        
        # 量子体积
        v_Q = lambda_th**3
        
        # 熵
        S = self.N * self.k_B * (np.log(self.V / (self.N * v_Q)) + 2.5)
        
        return S
    
    def heat_capacity_V(self) -> float:
        """
        定容热容
        
        C_V = (3/2) Nk_B
        """
        return 1.5 * self.N * self.k_B
    
    def heat_capacity_P(self) -> float:
        """
        定压热容
        
        C_P = (5/2) Nk_B = C_V + Nk_B
        """
        return 2.5 * self.N * self.k_B
    
    def adiabatic_index(self) -> float:
        """
        绝热指数
        
        γ = C_P / C_V
        """
        return self.heat_capacity_P() / self.heat_capacity_V()
    
    def maxwell_boltzmann_distribution(self, 
                                       velocities: np.ndarray) -> np.ndarray:
        """
        Maxwell-Boltzmann 速度分布
        
        f(v) = (m/2πk_BT)^(3/2) exp(-mv²/2k_BT)
        """
        factor = (self.m / (2 * np.pi * self.k_B * self.T))**(3/2)
        exponent = -self.m * velocities**2 / (2 * self.k_B * self.T)
        
        return factor * np.exp(exponent)
    
    def mean_speed(self) -> float:
        """
        平均速率
        
        ⟨v⟩ = √(8k_BT/πm)
        """
        return np.sqrt(8 * self.k_B * self.T / (np.pi * self.m))
    
    def rms_speed(self) -> float:
        """
        方均根速率
        
        v_rms = √(3k_BT/m)
        """
        return np.sqrt(3 * self.k_B * self.T / self.m)
    
    def most_probable_speed(self) -> float:
        """
        最概然速率
        
        v_p = √(2k_BT/m)
        """
        return np.sqrt(2 * self.k_B * self.T / self.m)


class VanDerWaalsGas(IdealGas):
    """
    van der Waals 气体
    
    考虑分子间相互作用和分子体积
    
    状态方程: (P + a/V²)(V - b) = Nk_B T
    
    体现 ±0: 吸引力(负) vs 排斥力(正)
    """
    
    def __init__(self,
                 n_particles: int,
                 temperature: float,
                 volume: float,
                 a: float,
                 b: float):
        """
        Args:
            a: 吸引力参数
            b: 排斥体积参数
        """
        super().__init__(n_particles, temperature, volume)
        self.a = a
        self.b = b
    
    def pressure(self) -> float:
        """
        van der Waals 压强
        
        P = Nk_B T/(V - Nb) - a(N/V)²
        """
        if self.V <= self.N * self.b:
            return float('inf')  # 物理不允许
        
        P_kinetic = self.N * self.k_B * self.T / (self.V - self.N * self.b)
        P_attractive = self.a * (self.N / self.V)**2
        
        return P_kinetic - P_attractive
    
    def critical_point(self) -> Tuple[float, float, float]:
        """
        临界点 (T_c, V_c, P_c)
        
        T_c = 8a/(27k_B b)
        V_c = 3Nb
        P_c = a/(27b²)
        """
        T_c = 8 * self.a / (27 * self.k_B * self.b)
        V_c = 3 * self.N * self.b
        P_c = self.a / (27 * self.b**2)
        
        return (T_c, V_c, P_c)
    
    def compressibility_factor(self) -> float:
        """
        压缩因子
        
        Z = PV/(Nk_BT)
        
        理想气体: Z = 1
        """
        P = self.pressure()
        return P * self.V / (self.N * self.k_B * self.T)


# ============================================================================
# Level 3: 量子统计
# ============================================================================

class QuantumStatistics(ABC):
    """
    量子统计抽象基类
    
    全同粒子的两种统计：
    - 玻色-Einstein 统计（整数自旋）
    - 费米-Dirac 统计（半整数自旋）
    """
    
    @abstractmethod
    def occupation_number(self, energy: float, 
                         temperature: float, 
                         chemical_potential: float) -> float:
        """单态平均占据数"""
        pass


class BoseEinsteinStatistics(QuantumStatistics):
    """
    玻色-Einstein 统计
    
    玻色子：可以多个粒子占据同一态
    
    例：光子、He-4、介子
    """
    
    def occupation_number(self, 
                         energy: float,
                         temperature: float,
                         chemical_potential: float = 0.0) -> float:
        """
        玻色-Einstein 分布
        
        n_BE(ε) = 1 / (e^((ε-μ)/k_BT) - 1)
        
        对于光子：μ = 0
        """
        k_B = PhysicalConstants.K_B_REDUCED
        beta = 1.0 / (k_B * temperature)
        
        exponent = beta * (energy - chemical_potential)
        
        if exponent > 100:  # 避免溢出
            return 0.0
        
        denominator = np.exp(exponent) - 1
        
        if abs(denominator) < 1e-10:
            return float('inf')  # Bose-Einstein 凝聚
        
        return 1.0 / denominator
    
    def bose_einstein_condensation_temperature(self,
                                              density: float,
                                              mass: float) -> float:
        """
        BEC 临界温度
        
        T_c = (2πℏ²/mk_B)(n/ζ(3/2))^(2/3)
        
        其中 ζ(3/2) ≈ 2.612 是 Riemann zeta 函数
        """
        h = PhysicalConstants.H
        k_B = PhysicalConstants.K_B_REDUCED
        
        zeta_3_2 = 2.612
        
        T_c = (2 * np.pi * (h / (2 * np.pi))**2 / (mass * k_B)) * \
              (density / zeta_3_2)**(2/3)
        
        return T_c


class FermiDiracStatistics(QuantumStatistics):
    """
    费米-Dirac 统计
    
    费米子：泡利不相容原理，每态最多一个粒子
    
    例：电子、质子、中子、He-3
    """
    
    def occupation_number(self,
                         energy: float,
                         temperature: float,
                         chemical_potential: float) -> float:
        """
        费米-Dirac 分布
        
        n_FD(ε) = 1 / (e^((ε-μ)/k_BT) + 1)
        """
        k_B = PhysicalConstants.K_B_REDUCED
        beta = 1.0 / (k_B * temperature)
        
        exponent = beta * (energy - chemical_potential)
        
        if exponent > 100:
            return 0.0
        elif exponent < -100:
            return 1.0
        
        return 1.0 / (np.exp(exponent) + 1)
    
    def fermi_energy(self, density: float, mass: float) -> float:
        """
        Fermi 能量（T=0 时的化学势）
        
        ε_F = (ℏ²/2m)(3π²n)^(2/3)
        """
        hbar = PhysicalConstants.HBAR
        
        return (hbar**2 / (2 * mass)) * (3 * np.pi**2 * density)**(2/3)
    
    def fermi_temperature(self, density: float, mass: float) -> float:
        """
        Fermi 温度
        
        T_F = ε_F / k_B
        """
        k_B = PhysicalConstants.K_B_REDUCED
        E_F = self.fermi_energy(density, mass)
        
        return E_F / k_B
    
    def heat_capacity_low_temperature(self,
                                     temperature: float,
                                     density: float,
                                     mass: float) -> float:
        """
        低温热容（T << T_F）
        
        C_V = (π²/2) n k_B (T/T_F)
        """
        k_B = PhysicalConstants.K_B_REDUCED
        T_F = self.fermi_temperature(density, mass)
        
        return (np.pi**2 / 2) * density * k_B * (temperature / T_F)


class BlackBodyRadiation:
    """
    黑体辐射
    
    光子气体（玻色子，μ=0）
    
    体现 ±0: 发射 ↔ 吸收 的平衡
    """
    
    def __init__(self, temperature: float):
        """
        Args:
            temperature: 黑体温度 T
        """
        self.T = temperature
        self.k_B = PhysicalConstants.K_B
        self.h = PhysicalConstants.H
        self.c = PhysicalConstants.C
    
    def planck_distribution(self, frequency: float) -> float:
        """
        Planck 谱分布
        
        u(ν) = (8πh/c³) ν³ / (e^(hν/k_BT) - 1)
        
        能量密度（单位频率）
        """
        if frequency <= 0:
            return 0.0
        
        factor = (8 * np.pi * self.h / self.c**3) * frequency**3
        exponent = self.h * frequency / (self.k_B * self.T)
        
        if exponent > 100:
            return 0.0
        
        denominator = np.exp(exponent) - 1
        
        if denominator < 1e-10:
            return 0.0
        
        return factor / denominator
    
    def wien_displacement_law(self) -> float:
        """
        Wien 位移定律
        
        λ_max T = 2.898 × 10⁻³ m·K
        """
        return 2.898e-3 / self.T
    
    def stefan_boltzmann_law(self) -> float:
        """
        Stefan-Boltzmann 定律
        
        j = σT⁴
        
        总辐射功率密度
        """
        sigma = PhysicalConstants.SIGMA
        return sigma * self.T**4
    
    def total_energy_density(self) -> float:
        """
        总能量密度
        
        u = (π²/15)(k_B T)⁴/(ℏc)³
        """
        hbar = PhysicalConstants.HBAR
        
        return (np.pi**2 / 15) * (self.k_B * self.T)**4 / (hbar * self.c)**3
    
    def photon_number_density(self) -> float:
        """
        光子数密度
        
        n = (2ζ(3)/π²)(k_B T/ℏc)³
        """
        hbar = PhysicalConstants.HBAR
        zeta_3 = 1.202  # Riemann zeta(3)
        
        return (2 * zeta_3 / np.pi**2) * (self.k_B * self.T / (hbar * self.c))**3
    
    def entropy_density(self) -> float:
        """
        熵密度
        
        s = (4/3) u/T
        """
        u = self.total_energy_density()
        return (4/3) * u / self.T


# ============================================================================
# Level 4: 热力学定律
# ============================================================================

class ThermodynamicLaws:
    """
    热力学四大定律
    
    基于 ±0 的理解:
    - 第零定律: 热平衡（传递性）
    - 第一定律: 能量守恒（正负均衡）
    - 第二定律: 熵增（趋向最大均衡）
    - 第三定律: 绝对零度不可达（基态唯一性）
    """
    
    @staticmethod
    def first_law(internal_energy_change: float,
                  heat: float,
                  work: float) -> bool:
        """
        第一定律
        
        ΔU = Q - W
        
        验证能量守恒
        """
        return abs(internal_energy_change - (heat - work)) < 1e-10
    
    @staticmethod
    def second_law_clausius(heat_in: float,
                           heat_out: float,
                           T_hot: float,
                           T_cold: float) -> bool:
        """
        第二定律（Clausius 表述）
        
        热量不能自发从冷物体流向热物体
        
        ΔS_total = Q_out/T_cold - Q_in/T_hot ≥ 0
        """
        if T_hot <= T_cold:
            return False
        
        dS_hot = -heat_in / T_hot
        dS_cold = heat_out / T_cold
        dS_total = dS_hot + dS_cold
        
        return dS_total >= -1e-10
    
    @staticmethod
    def second_law_kelvin(work: float, heat_absorbed: float) -> bool:
        """
        第二定律（Kelvin-Planck 表述）
        
        不可能从单一热源吸热全部转化为功
        
        效率 η < 1
        """
        if heat_absorbed <= 0:
            return True
        
        efficiency = work / heat_absorbed
        return efficiency < 1.0 + 1e-10
    
    @staticmethod
    def carnot_efficiency(T_hot: float, T_cold: float) -> float:
        """
        Carnot 效率（可逆热机的最高效率）
        
        η_C = 1 - T_cold/T_hot
        
        体现 ±0: 温差越大，可用能量越多
        """
        if T_hot <= T_cold:
            return 0.0
        
        return 1.0 - T_cold / T_hot
    
    @staticmethod
    def third_law_entropy(temperature: float, 
                         S_0: float = 0.0) -> float:
        """
        第三定律
        
        S(T→0) → S_0 (通常 S_0 = 0)
        
        绝对零度时，完美晶体的熵为零
        """
        if temperature < 1e-10:
            return S_0
        
        # 简化：S ∝ T³ (Debye 模型)
        return S_0 + temperature**3


class HeatEngine:
    """
    热机
    
    循环过程：吸热 → 做功 → 放热
    
    体现 ±0: 能量转换中的正负流动
    """
    
    def __init__(self, T_hot: float, T_cold: float):
        """
        Args:
            T_hot: 高温热源温度
            T_cold: 低温热源温度
        """
        self.T_hot = T_hot
        self.T_cold = T_cold
    
    def carnot_cycle(self, Q_in: float) -> Tuple[float, float, float]:
        """
        Carnot 循环
        
        Returns:
            (做功 W, 放热 Q_out, 效率 η)
        """
        eta = ThermodynamicLaws.carnot_efficiency(self.T_hot, self.T_cold)
        W = eta * Q_in
        Q_out = Q_in - W
        
        return (W, Q_out, eta)
    
    def otto_cycle(self, 
                   compression_ratio: float,
                   gamma: float = 1.4) -> float:
        """
        Otto 循环（汽油机）
        
        η = 1 - 1/r^(γ-1)
        
        Args:
            compression_ratio: 压缩比 r
            gamma: 绝热指数
        """
        return 1.0 - 1.0 / (compression_ratio**(gamma - 1))
    
    def diesel_cycle(self,
                    compression_ratio: float,
                    cutoff_ratio: float,
                    gamma: float = 1.4) -> float:
        """
        Diesel 循环（柴油机）
        
        η = 1 - (1/r^(γ-1)) × (r_c^γ - 1)/(γ(r_c - 1))
        """
        r = compression_ratio
        r_c = cutoff_ratio
        
        term1 = 1.0 / (r**(gamma - 1))
        term2 = (r_c**gamma - 1) / (gamma * (r_c - 1))
        
        return 1.0 - term1 * term2


class Refrigerator:
    """
    制冷机
    
    逆向热机：做功 → 从低温吸热 → 向高温放热
    """
    
    def __init__(self, T_hot: float, T_cold: float):
        self.T_hot = T_hot
        self.T_cold = T_cold
    
    def coefficient_of_performance(self) -> float:
        """
        制冷系数 (COP)
        
        COP = Q_cold / W = T_cold / (T_hot - T_cold)
        
        (理想 Carnot 制冷机)
        """
        if self.T_hot <= self.T_cold:
            return float('inf')
        
        return self.T_cold / (self.T_hot - self.T_cold)
    
    def heat_pump_performance(self) -> float:
        """
        热泵系数
        
        COP_HP = Q_hot / W = T_hot / (T_hot - T_cold)
        """
        if self.T_hot <= self.T_cold:
            return float('inf')
        
        return self.T_hot / (self.T_hot - self.T_cold)


# ============================================================================
# Level 5: 相变与临界现象
# ============================================================================

class PhaseTransition:
    """
    相变
    
    物质在不同相之间的转换
    
    基于 ±0 的理解:
    - 有序相 ↔ 无序相
    - 对称破缺
    - 临界点：两相均衡
    """
    
    def __init__(self, 
                 T_critical: float,
                 order_parameter_func: Callable[[float], float]):
        """
        Args:
            T_critical: 临界温度
            order_parameter_func: 序参量函数 m(T)
        """
        self.T_c = T_critical
        self.order_param = order_parameter_func
    
    def is_first_order(self, T: float) -> bool:
        """
        一级相变：序参量跳变
        
        例：水-冰，液-气
        """
        dT = 1e-6
        m1 = self.order_param(T - dT)
        m2 = self.order_param(T + dT)
        
        return abs(m2 - m1) > 1e-3
    
    def is_second_order(self, T: float) -> bool:
        """
        二级相变：序参量连续，导数跳变
        
        例：铁磁-顺磁转变
        """
        return not self.is_first_order(T)
    
    def critical_exponents(self) -> Dict[str, float]:
        """
        临界指数
        
        - α: 热容 C ∝ |T - T_c|^(-α)
        - β: 序参量 m ∝ |T - T_c|^β
        - γ: 磁化率 χ ∝ |T - T_c|^(-γ)
        - δ: 临界等温线 m ∝ h^(1/δ)
        - ν: 关联长度 ξ ∝ |T - T_c|^(-ν)
        """
        # 平均场理论值
        return {
            'alpha': 0.0,
            'beta': 0.5,
            'gamma': 1.0,
            'delta': 3.0,
            'nu': 0.5
        }


class IsingModel:
    """
    Ising 模型
    
    最简单的相互作用自旋系统
    
    H = -J Σ_{<i,j>} s_i s_j - h Σ_i s_i
    
    体现 ±0: 自旋向上(+1) ↔ 自旋向下(-1)
    """
    
    def __init__(self,
                 lattice_size: int,
                 coupling: float,
                 external_field: float = 0.0):
        """
        Args:
            lattice_size: 格点数（一维）
            coupling: 耦合常数 J
            external_field: 外场 h
        """
        self.L = lattice_size
        self.J = coupling
        self.h = external_field
        
        # 初始化自旋配置（随机）
        self.spins = np.random.choice([-1, 1], size=self.L)
    
    def energy(self) -> float:
        """
        系统能量
        
        E = -J Σ s_i s_{i+1} - h Σ s_i
        """
        # 近邻相互作用（周期边界条件）
        interaction = 0.0
        for i in range(self.L):
            j = (i + 1) % self.L
            interaction += self.spins[i] * self.spins[j]
        
        E_interaction = -self.J * interaction
        E_field = -self.h * np.sum(self.spins)
        
        return E_interaction + E_field
    
    def magnetization(self) -> float:
        """
        磁化强度（序参量）
        
        m = (1/N) Σ s_i
        """
        return np.mean(self.spins)
    
    def monte_carlo_step(self, temperature: float, k_B: float = 1.0):
        """
        Monte Carlo 更新（Metropolis 算法）
        """
        beta = 1.0 / (k_B * temperature)
        
        for _ in range(self.L):
            # 随机选择一个格点
            i = np.random.randint(self.L)
            
            # 计算翻转能量变化
            left = (i - 1) % self.L
            right = (i + 1) % self.L
            
            dE = 2 * self.spins[i] * (self.J * (self.spins[left] + self.spins[right]) + self.h)
            
            # Metropolis 接受准则
            if dE < 0 or np.random.rand() < np.exp(-beta * dE):
                self.spins[i] *= -1
    
    def simulate(self, 
                temperature: float,
                n_steps: int = 1000,
                thermalization: int = 100) -> Tuple[float, float]:
        """
        Monte Carlo 模拟
        
        Returns:
            (平均磁化强度, 平均能量)
        """
        # 热化
        for _ in range(thermalization):
            self.monte_carlo_step(temperature)
        
        # 测量
        magnetizations = []
        energies = []
        
        for _ in range(n_steps):
            self.monte_carlo_step(temperature)
            magnetizations.append(abs(self.magnetization()))
            energies.append(self.energy())
        
        return (np.mean(magnetizations), np.mean(energies))
    
    @staticmethod
    def critical_temperature_2d(J: float, k_B: float = 1.0) -> float:
        """
        2D Ising 模型临界温度（精确解）
        
        k_B T_c = (2/ln(1+√2)) J ≈ 2.269 J
        """
        return (2.0 / np.log(1 + np.sqrt(2))) * J / k_B


# ============================================================================
# Level 6: 非平衡态热力学
# ============================================================================

class NonEquilibriumThermodynamics:
    """
    非平衡态热力学
    
    系统偏离平衡态的行为
    
    基于 ±0: 从不均衡向均衡演化
    """
    
    @staticmethod
    def heat_conduction(
        temperature_gradient: np.ndarray,
        thermal_conductivity: float
    ) -> np.ndarray:
        """
        Fourier 定律
        
        j_q = -κ ∇T
        
        热流从高温流向低温
        """
        return -thermal_conductivity * temperature_gradient
    
    @staticmethod
    def diffusion(
        concentration_gradient: np.ndarray,
        diffusion_coefficient: float
    ) -> np.ndarray:
        """
        Fick 定律
        
        j_n = -D ∇n
        
        粒子从高浓度扩散到低浓度
        """
        return -diffusion_coefficient * concentration_gradient
    
    @staticmethod
    def viscous_flow(
        velocity_gradient: np.ndarray,
        viscosity: float
    ) -> np.ndarray:
        """
        Newton 粘性定律
        
        τ = η ∂v/∂x
        """
        return viscosity * velocity_gradient
    
    @staticmethod
    def onsager_reciprocal_relations(
        fluxes: np.ndarray,
        forces: np.ndarray,
        L_matrix: np.ndarray
    ) -> bool:
        """
        Onsager 倒易关系
        
        J_i = Σ_j L_ij X_j
        L_ij = L_ji (线性响应区)
        """
        # 验证对称性
        return np.allclose(L_matrix, L_matrix.T)
    
    @staticmethod
    def entropy_production_rate(
        fluxes: np.ndarray,
        forces: np.ndarray
    ) -> float:
        """
        熵产生率
        
        dS/dt = Σ_i J_i X_i ≥ 0
        
        体现不可逆性
        """
        return np.dot(fluxes, forces)


class BoltzmannEquation:
    """
    Boltzmann 方程
    
    描述气体分子速度分布的演化
    
    ∂f/∂t + v·∇f + (F/m)·∇_v f = (∂f/∂t)_coll
    """
    
    def __init__(self, 
                 collision_time: float,
                 mass: float):
        """
        Args:
            collision_time: 碰撞弛豫时间 τ
            mass: 粒子质量 m
        """
        self.tau = collision_time
        self.m = mass
    
    def relaxation_time_approximation(self,
                                     f: np.ndarray,
                                     f_eq: np.ndarray) -> np.ndarray:
        """
        弛豫时间近似
        
        (∂f/∂t)_coll = -(f - f_eq)/τ
        """
        return -(f - f_eq) / self.tau
    
    def h_theorem(self, 
                 f: np.ndarray,
                 velocities: np.ndarray) -> float:
        """
        H 定理
        
        H = ∫ f ln f d³v
        dH/dt ≤ 0
        
        证明熵增原理
        """
        # 避免 log(0)
        f_safe = np.where(f > 1e-10, f, 1e-10)
        
        H = np.sum(f_safe * np.log(f_safe))
        
        return H


# ============================================================================
# 综合演示系统
# ============================================================================

class ThermodynamicsDemo:
    """热力学综合演示"""
    
    @staticmethod
    def demonstrate_ensembles():
        """演示统计系综"""
        print("\n" + "=" * 70)
        print("统计系综")
        print("=" * 70)
        
        # 创建能级
        energies = [0, 1, 2, 3, 4]
        levels = [EnergyLevel(E, degeneracy=2*i+1) 
                 for i, E in enumerate(energies)]
        
        # 正则系综
        T = 1.0
        canonical = CanonicalEnsemble(T, volume=1.0, particle_number=10, 
                                     energy_levels=levels)
        
        Z = canonical.partition_function()
        E_avg = canonical.average_energy()
        S = canonical.entropy()
        F = canonical.free_energy()
        C = canonical.heat_capacity()
        
        print(f"\n正则系综 (T={T}):")
        print(f"  配分函数 Z = {Z:.4f}")
        print(f"  平均能量 ⟨E⟩ = {E_avg:.4f}")
        print(f"  熵 S = {S:.4f}")
        print(f"  自由能 F = {F:.4f}")
        print(f"  热容 C = {C:.4f}")
        
        # 概率分布
        probs = canonical.probability_distribution()
        print(f"\n能级占据概率:")
        for i, (level, prob) in enumerate(zip(levels, probs)):
            print(f"  E_{i} = {level.energy}: P = {prob:.4f}")
    
    @staticmethod
    def demonstrate_ideal_gas():
        """演示理想气体"""
        print("\n" + "=" * 70)
        print("理想气体")
        print("=" * 70)
        
        N = 1000
        T = 300  # K
        V = 1.0  # m³
        
        gas = IdealGas(N, T, V)
        
        P = gas.pressure()
        U = gas.internal_energy()
        S = gas.entropy()
        
        print(f"\n理想气体 (N={N}, T={T}K, V={V}m³):")
        print(f"  压强 P = {P:.2e} Pa")
        print(f"  内能 U = {U:.2e} J")
        print(f"  熵 S = {S:.2e} J/K")
        
        # 速度分布
        v_mean = gas.mean_speed()
        v_rms = gas.rms_speed()
        v_p = gas.most_probable_speed()
        
        print(f"\n速度统计:")
        print(f"  平均速率 ⟨v⟩ = {v_mean:.2f} m/s")
        print(f"  方均根速率 v_rms = {v_rms:.2f} m/s")
        print(f"  最概然速率 v_p = {v_p:.2f} m/s")
    
    @staticmethod
    def demonstrate_quantum_statistics():
        """演示量子统计"""
        print("\n" + "=" * 70)
        print("量子统计")
        print("=" * 70)
        
        # 费米-Dirac
        fermi = FermiDiracStatistics()
        
        density = 1e28  # m⁻³
        mass = 9.109e-31  # 电子质量
        
        E_F = fermi.fermi_energy(density, mass)
        T_F = fermi.fermi_temperature(density, mass)
        
        print(f"\n费米气体:")
        print(f"  密度 n = {density:.2e} m⁻³")
        print(f"  Fermi 能量 E_F = {E_F:.2e} J")
        print(f"  Fermi 温度 T_F = {T_F:.2e} K")
        
        # 不同温度下的占据数
        energies = np.linspace(0, 2*E_F, 100)
        for T in [0.1*T_F, T_F, 10*T_F]:
            print(f"\n  T = {T:.2e} K:")
            n_E_F = fermi.occupation_number(E_F, T, E_F)
            print(f"    n(E_F) = {n_E_F:.4f}")
        
        # 玻色-Einstein
        bose = BoseEinsteinStatistics()
        
        T_BEC = bose.bose_einstein_condensation_temperature(density, mass)
        print(f"\nBose 气体:")
        print(f"  BEC 临界温度 T_c = {T_BEC:.2e} K")
    
    @staticmethod
    def demonstrate_blackbody():
        """演示黑体辐射"""
        print("\n" + "=" * 70)
        print("黑体辐射")
        print("=" * 70)
        
        T = 5778  # K (太阳表面温度)
        blackbody = BlackBodyRadiation(T)
        
        lambda_max = blackbody.wien_displacement_law()
        j_total = blackbody.stefan_boltzmann_law()
        
        print(f"\n黑体 (T={T}K):")
        print(f"  Wien 位移定律 λ_max = {lambda_max*1e9:.1f} nm")
        print(f"  Stefan-Boltzmann J = {j_total:.2e} W/m²")
        
        # Planck 分布
        print(f"\nPlanck 谱分布:")
        wavelengths = [400e-9, 550e-9, 700e-9]  # 可见光
        for wl in wavelengths:
            nu = blackbody.c / wl
            u_nu = blackbody.planck_distribution(nu)
            print(f"  λ = {wl*1e9:.0f} nm: u = {u_nu:.2e}")
    
    @staticmethod
    def demonstrate_heat_engine():
        """演示热机"""
        print("\n" + "=" * 70)
        print("热机")
        print("=" * 70)
        
        T_hot = 600  # K
        T_cold = 300  # K
        
        engine = HeatEngine(T_hot, T_cold)
        
        # Carnot 循环
        Q_in = 1000  # J
        W, Q_out, eta = engine.carnot_cycle(Q_in)
        
        print(f"\nCarnot 热机:")
        print(f"  高温热源 T_H = {T_hot} K")
        print(f"  低温热源 T_C = {T_cold} K")
        print(f"  吸热 Q_in = {Q_in} J")
        print(f"  做功 W = {W:.2f} J")
        print(f"  放热 Q_out = {Q_out:.2f} J")
        print(f"  效率 η = {eta:.2%}")
        
        # Otto 循环
        eta_otto = engine.otto_cycle(compression_ratio=9)
        print(f"\nOtto 循环 (压缩比=9):")
        print(f"  效率 η = {eta_otto:.2%}")
        
        # 制冷机
        refrigerator = Refrigerator(T_hot, T_cold)
        COP = refrigerator.coefficient_of_performance()
        
        print(f"\n制冷机:")
        print(f"  制冷系数 COP = {COP:.2f}")
    
    @staticmethod
    def demonstrate_phase_transition():
        """演示相变"""
        print("\n" + "=" * 70)
        print("相变与临界现象")
        print("=" * 70)
        
        # Ising 模型
        L = 50
        J = 1.0
        ising = IsingModel(L, J)
        
        T_c = IsingModel.critical_temperature_2d(J)
        print(f"\n1D Ising 模型 (L={L}):")
        print(f"  临界温度 T_c = {T_c:.3f}")
        
        # 不同温度下的磁化强度
        temperatures = [0.5*T_c, T_c, 2*T_c]
        print(f"\nMonte Carlo 模拟:")
        
        for T in temperatures:
            m_avg, E_avg = ising.simulate(T, n_steps=1000)
            print(f"  T/T_c = {T/T_c:.2f}: m = {m_avg:.4f}, E = {E_avg:.2f}")


def main():
    """主函数：完整演示"""
    
    print("=" * 70)
    print("DaoMath Thermodynamics & Statistical Physics")
    print("从 ±0 构造热力学")
    print("=" * 70)
    
    # 1. 统计系综
    ThermodynamicsDemo.demonstrate_ensembles()
    
    # 2. 理想气体
    ThermodynamicsDemo.demonstrate_ideal_gas()
    
    # 3. 量子统计
    ThermodynamicsDemo.demonstrate_quantum_statistics()
    
    # 4. 黑体辐射
    ThermodynamicsDemo.demonstrate_blackbody()
    
    # 5. 热机
    ThermodynamicsDemo.demonstrate_heat_engine()
    
    # 6. 相变
    ThermodynamicsDemo.demonstrate_phase_transition()
    
    print("\n" + "=" * 70)
    print("热力学理论构造完成！")
    print("=" * 70)


if __name__ == "__main__":
    main()
